# Load libraries
library(data.table)
library(devtools)
library(presto)
library(glmGamPoi)
library(sctransform)
library(Seurat)
library(tidyverse)
library(miQC)
library(SeuratWrappers)
library(flexmix)
library(SingleCellExperiment)
library(SummarizedExperiment)
library(readxl)
library(Matrix)
library(speckle)
library(scater)
library(patchwork)
library(vctrs)
library(harmony)
library(scDblFinder)
library(ggridges)
library(scGate)
library(ggsankey)
library(ggpubr)
# Set global options for Seurat v5 objects
options(Seurat.object.assay.version = 'v5')Roider DLBCL scRNA-seq reprocessing
Rendered 12/12/2025
Set up Seurat workspace
Import 10x data and create Seurat objects with sample names prepended to cell barcodes
DLBCL1 <- Read10X("DLBCL1",strip.suffix = T)
DLBCL2 <- Read10X("DLBCL2",strip.suffix = T)
DLBCL3 <- Read10X("DLBCL3",strip.suffix = T)
FL1 <- Read10X("FL1",strip.suffix = T)
FL2 <- Read10X("FL2",strip.suffix = T)
FL3 <- Read10X("FL3",strip.suffix = T)
FL4 <- Read10X("FL4",strip.suffix = T)
rLN1 <- Read10X("rLN1",strip.suffix = T)
rLN2 <- Read10X("rLN2",strip.suffix = T)
rLN3 <- Read10X("rLN3",strip.suffix = T)
tFL1 <- Read10X("tFL1",strip.suffix = T)
tFL2 <- Read10X("tFL2",strip.suffix = T)
DLBCL1.seurat <- CreateSeuratObject(counts = DLBCL1)Warning: Feature names cannot have underscores ('_'), replacing with dashes
('-')DLBCL1.seurat <- RenameCells(DLBCL1.seurat,add.cell.id = "DLBCL1")
DLBCL2.seurat <- CreateSeuratObject(counts = DLBCL2)Warning: Feature names cannot have underscores ('_'), replacing with dashes
('-')DLBCL2.seurat <- RenameCells(DLBCL2.seurat,add.cell.id = "DLBCL2")
DLBCL3.seurat <- CreateSeuratObject(counts = DLBCL3)
DLBCL3.seurat <- RenameCells(DLBCL3.seurat,add.cell.id = "DLBCL3")
FL1.seurat <- CreateSeuratObject(counts = FL1)Warning: Feature names cannot have underscores ('_'), replacing with dashes
('-')FL1.seurat <- RenameCells(FL1.seurat,add.cell.id = "FL1")
FL2.seurat <- CreateSeuratObject(counts = FL2)Warning: Feature names cannot have underscores ('_'), replacing with dashes
('-')FL2.seurat <- RenameCells(FL2.seurat,add.cell.id = "FL2")
FL3.seurat <- CreateSeuratObject(counts = FL3)Warning: Feature names cannot have underscores ('_'), replacing with dashes
('-')FL3.seurat <- RenameCells(FL3.seurat,add.cell.id = "FL3")
FL4.seurat <- CreateSeuratObject(counts = FL4)
FL4.seurat <- RenameCells(FL4.seurat,add.cell.id = "FL4")
rLN1.seurat <- CreateSeuratObject(counts = rLN1)Warning: Feature names cannot have underscores ('_'), replacing with dashes
('-')rLN1.seurat <- RenameCells(rLN1.seurat,add.cell.id = "rLN1")
rLN2.seurat <- CreateSeuratObject(counts = rLN2)
rLN2.seurat <- RenameCells(rLN2.seurat,add.cell.id = "rLN2")
rLN3.seurat <- CreateSeuratObject(counts = rLN3)
rLN3.seurat <- RenameCells(rLN3.seurat,add.cell.id = "rLN3")
tFL1.seurat <- CreateSeuratObject(counts = tFL1)Warning: Feature names cannot have underscores ('_'), replacing with dashes
('-')tFL1.seurat <- RenameCells(tFL1.seurat,add.cell.id = "tFL1")
tFL2.seurat <- CreateSeuratObject(counts = tFL2)
tFL2.seurat <- RenameCells(tFL2.seurat,add.cell.id = "tFL2")Retrieve published cell-cluster metadata
This is seemingly on a downsampled set of cells, perhaps done so for easier plotting in the shiny app
download.file("https://github.com/anders-biostat/single-cell-lymphoma-explorer/raw/master/metaInfo.rds",
destfile = "Roider_DLBCL_metaInfo.rds")
roider.meta <- readRDS("Roider_DLBCL_metaInfo.rds")
head(roider.meta$bcells) UMAP1 UMAP2 Subset
tFL2_GATTCAGAGCTGGAAC -3.726201 -7.181145 5
tFL2_AACACGTCACACAGAG -3.260301 -5.485855 5
tFL2_ACATCAGGTTCCGGCA -4.213102 -6.762125 5
tFL2_CCCAGTTTCACCGGGT -3.345386 -5.145964 5
tFL2_ACTTGTTTCGTGGACC -4.179855 -6.458187 5
tFL2_CTCGTCATCACCGGGT -2.707596 -6.752877 5Merge objects
merged.roider <- merge(x = DLBCL1.seurat, y=c(DLBCL2.seurat, DLBCL3.seurat, FL1.seurat, FL2.seurat, FL3.seurat, FL4.seurat, rLN1.seurat, rLN2.seurat, rLN3.seurat, tFL1.seurat, tFL2.seurat))
dim(merged.roider)[1] 45068 41786QC filter
merged.roider <- subset(merged.roider, subset = nFeature_RNA > 250)
dim(merged.roider)[1] 45068 41645merged.roider <- PercentageFeatureSet(merged.roider, pattern = "^MT-", col.name = "percent.mt")
merged.roider <- RunMiQC(merged.roider,
percent.mt = "percent.mt",
nFeature_RNA = "nFeature_RNA",
posterior.cutoff = 0.7,
model.slot = "flexmix_model")Warning: Adding a command log without an assay associated with itmerged.roider <- subset(merged.roider, miQC.keep == "keep")
dim(merged.roider)[1] 45068 40073data.frame(table(str_replace_all(colnames(merged.roider),"_.+",""))) Var1 Freq
1 DLBCL1 4008
2 DLBCL2 4927
3 DLBCL3 2412
4 FL1 3668
5 FL2 5090
6 FL3 4686
7 FL4 3147
8 rLN1 3082
9 rLN2 2028
10 rLN3 2779
11 tFL1 1998
12 tFL2 2248Add meta data
disease <- str_replace_all(colnames(merged.roider),"\\d_.+","")
sample <- str_replace_all(colnames(merged.roider),"_.+","")
merged.roider <- AddMetaData(merged.roider,disease,col.name="Disease")
merged.roider <- AddMetaData(merged.roider,sample,col.name="Sample")Join then re-split RNA counts layers by Sample
merged.roider[['RNA']] <- JoinLayers(merged.roider[['RNA']])
merged.roider[["RNA"]] <- split(merged.roider[["RNA"]], f = merged.roider$Sample)Normalize and scale data
Regress out mitochondrial contribution
merged.roider <- NormalizeData(merged.roider)
merged.roider <- FindVariableFeatures(merged.roider,
assay="RNA",
layer="data",
selection.method = "vst",
nfeatures = 5000)
merged.roider <- ScaleData(merged.roider, vars.to.regress = "percent.mt")Run PCA
merged.roider <- RunPCA(merged.roider, npcs = 200, verbose = FALSE)
ElbowPlot(merged.roider, ndims = 200, reduction = "pca")
Plot unintegrated UMAP
merged.roider <- RunUMAP(merged.roider,
reduction = "pca",
dims = 1:30,
reduction.name = "umap.unintegrated")Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session10:55:58 UMAP embedding parameters a = 0.9922 b = 1.112Found more than one class "dist" in cache; using the first, from namespace 'BiocGenerics'Also defined by 'spam'10:55:58 Read 40073 rows and found 30 numeric columns10:55:58 Using Annoy for neighbor search, n_neighbors = 30Found more than one class "dist" in cache; using the first, from namespace 'BiocGenerics'Also defined by 'spam'10:55:58 Building Annoy index with metric = cosine, n_trees = 500% 10 20 30 40 50 60 70 80 90 100%[----|----|----|----|----|----|----|----|----|----|**************************************************|
10:56:09 Writing NN index file to temp file /tmp/RtmpGN9GPk/file150e28686366a7
10:56:09 Searching Annoy index using 1 thread, search_k = 3000
10:56:25 Annoy recall = 100%
10:56:26 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
10:56:30 Initializing from normalized Laplacian + noise (using RSpectra)
10:56:41 Commencing optimization for 200 epochs, with 1735844 positive edges
10:57:12 Optimization finishedDimPlot(merged.roider, reduction = "umap.unintegrated", group.by = c("Disease", "Sample"))
Call preliminary clusters for the purposes of doublet discrimination
merged.roider <- FindNeighbors(merged.roider, dims = 1:30, reduction = "pca")Computing nearest neighbor graphComputing SNNmerged.roider <- FindClusters(merged.roider,
resolution = 0.2,
algorithm = 2)Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 40073
Number of edges: 1491903
Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.9798
Number of communities: 21
Elapsed time: 10 secondsDimPlot(merged.roider, reduction = "umap.unintegrated", label = T)
Identify and remove doublets
This uses raw counts as input
Convert seurat to sce and check colData
# Join RNA layers
merged.roider[['RNA']] <- JoinLayers(merged.roider[['RNA']])
merged.roider.sce <- as.SingleCellExperiment(merged.roider, assay = "RNA")
merged.roider.sceclass: SingleCellExperiment
dim: 45068 40073
metadata(0):
assays(2): counts logcounts
rownames(45068): RP11-34P13.3 FAM138A ... AC145212.1 MAFIP
rowData names(0):
colnames(40073): DLBCL1_AAACCTGAGCCATCGC DLBCL1_AAACCTGAGCCTTGAT ...
tFL2_TTTGTCATCAGGCGAA tFL2_TTTGTCATCAGGTAAA
colData names(11): orig.ident nCount_RNA ... seurat_clusters ident
reducedDimNames(2): PCA UMAP.UNINTEGRATED
mainExpName: RNA
altExpNames(0):colData(merged.roider.sce)DataFrame with 40073 rows and 11 columns
orig.ident nCount_RNA nFeature_RNA percent.mt
<character> <numeric> <integer> <numeric>
DLBCL1_AAACCTGAGCCATCGC SeuratProject 9180 2007 4.65142
DLBCL1_AAACCTGAGCCTTGAT SeuratProject 16902 3280 1.36670
DLBCL1_AAACCTGAGTTCGCAT SeuratProject 4800 1490 5.06250
DLBCL1_AAACCTGCACCACCAG SeuratProject 17052 4003 2.18743
DLBCL1_AAACCTGGTTAAGATG SeuratProject 8116 2297 1.88517
... ... ... ... ...
tFL2_TTTGTCAGTCATCGGC SeuratProject 7075 2077 2.388693
tFL2_TTTGTCAGTTCACGGC SeuratProject 5901 1677 1.745467
tFL2_TTTGTCATCAACGCTA SeuratProject 6722 2240 2.424874
tFL2_TTTGTCATCAGGCGAA SeuratProject 4281 1714 0.981079
tFL2_TTTGTCATCAGGTAAA SeuratProject 5842 1760 1.763095
miQC.probability miQC.keep Disease Sample
<numeric> <character> <character> <character>
DLBCL1_AAACCTGAGCCATCGC 0.05811154 keep DLBCL DLBCL1
DLBCL1_AAACCTGAGCCTTGAT 0.01431377 keep DLBCL DLBCL1
DLBCL1_AAACCTGAGTTCGCAT 0.16824566 keep DLBCL DLBCL1
DLBCL1_AAACCTGCACCACCAG 0.00701990 keep DLBCL DLBCL1
DLBCL1_AAACCTGGTTAAGATG 0.00628381 keep DLBCL DLBCL1
... ... ... ... ...
tFL2_TTTGTCAGTCATCGGC 0.00509071 keep tFL tFL2
tFL2_TTTGTCAGTTCACGGC 0.00515423 keep tFL tFL2
tFL2_TTTGTCATCAACGCTA 0.00530276 keep tFL tFL2
tFL2_TTTGTCATCAGGCGAA 0.01042485 keep tFL tFL2
tFL2_TTTGTCATCAGGTAAA 0.00532506 keep tFL tFL2
RNA_snn_res.0.2 seurat_clusters ident
<factor> <factor> <factor>
DLBCL1_AAACCTGAGCCATCGC 1 1 1
DLBCL1_AAACCTGAGCCTTGAT 1 1 1
DLBCL1_AAACCTGAGTTCGCAT 1 1 1
DLBCL1_AAACCTGCACCACCAG 1 1 1
DLBCL1_AAACCTGGTTAAGATG 1 1 1
... ... ... ...
tFL2_TTTGTCAGTCATCGGC 10 10 10
tFL2_TTTGTCAGTTCACGGC 10 10 10
tFL2_TTTGTCATCAACGCTA 10 10 10
tFL2_TTTGTCATCAGGCGAA 5 5 5
tFL2_TTTGTCATCAGGTAAA 10 10 10Run scDblFinder
Set the dbr.sd very high to better allow thresholds to be set based on misclassification rates per sample
merged.roider.sce <- scDblFinder(merged.roider.sce,
samples = "Sample",
dbr.sd = 1,
clusters = "seurat_clusters")Inspect results
# Look at the classes
table(merged.roider.sce$seurat_clusters, merged.roider.sce$scDblFinder.class)
singlet doublet
0 6291 215
1 3932 62
2 3945 21
3 3366 82
4 3371 36
5 2866 75
6 2609 16
7 2561 9
8 2243 38
9 2163 13
10 1635 33
11 767 76
12 689 151
13 631 205
14 740 62
15 272 217
16 97 110
17 140 39
18 0 120
19 88 25
20 39 23table(merged.roider.sce$Sample, merged.roider.sce$scDblFinder.class)
singlet doublet
DLBCL1 3932 76
DLBCL2 4662 265
DLBCL3 2354 58
FL1 3532 136
FL2 4951 139
FL3 4434 252
FL4 3039 108
rLN1 2833 249
rLN2 1957 71
rLN3 2710 69
tFL1 1874 124
tFL2 2167 81# Look at the scores
summary(merged.roider.sce$scDblFinder.score) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000757 0.0020391 0.0033650 0.0572434 0.0138431 0.9999380 Save doublet classifications into main Seurat object
merged.roider$doublet_classification <- merged.roider.sce$scDblFinder.classCount singlets and doublets
table(merged.roider$doublet_classification)
singlet doublet
38445 1628 table(merged.roider$doublet_classification, merged.roider$seurat_clusters)
0 1 2 3 4 5 6 7 8 9 10 11 12 13
singlet 6291 3932 3945 3366 3371 2866 2609 2561 2243 2163 1635 767 689 631
doublet 215 62 21 82 36 75 16 9 38 13 33 76 151 205
14 15 16 17 18 19 20
singlet 740 272 97 140 0 88 39
doublet 62 217 110 39 120 25 23Plot singlets/doublets in UMAP space
DimPlot(merged.roider,reduction = "umap.unintegrated", group.by = "doublet_classification")
Subset object to remove doublets and count remaining cells
merged.roider.singlets <- subset(merged.roider, doublet_classification %in% c("singlet"))
dim(merged.roider.singlets)[1] 45068 38445# Count remaining cells per initial cluster
table(merged.roider.singlets$seurat_clusters)
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
6291 3932 3945 3366 3371 2866 2609 2561 2243 2163 1635 767 689 631 740 272
16 17 18 19 20
97 140 0 88 39 Set up gating strategy to focus only on B-cell lineage
We gate to remove: * CD3+ (T cells) * DERL3+/ITM2C+ (Plasma cells)
# Plot genes to be used in gating model
FeaturePlot(merged.roider.singlets,
features=c("CD79A", "CD79B",
"MS4A1", "BANK1",
"CD3E", "CD3D",
"CD2", "GZMK",
"DERL3", "ITM2C"),
order = T,
ncol = 5)
Set up gating model and run scGate
scGate_model <- gating_model(name = "Tcells",
signature = c("CD3E", "CD3D",
"CD2", "GZMK"),
negative = T,
positive = F,
level = 1
)
scGate_model <- gating_model(model = scGate_model,
name = "Plasma",
signature = c("DERL3", "ITM2C"),
negative = T,
positive = F,
level = 2
)
scGate_model <- gating_model(model = scGate_model,
name = "Bcells",
signature = c("CD79A", "CD79B",
"MS4A1", "BANK1"),
positive = T,
negative = F,
level = 3
)
merged.roider.singlets <- scGate(data = merged.roider.singlets,
model = scGate_model,
pos.thr = 0.25,
neg.thr = 0.15,
verbose = T,
ncores = 4,
pca.dim = 40,
reduction = "umap.unintegrated")Computing UCell scores for all signatures using RNA assay...# Visualize "Pure" and "Impure" cells
DimPlot(merged.roider.singlets, group.by = "is.pure")
# Subset and remove non-B cells
merged.roider.singlets <- subset(merged.roider.singlets, is.pure %in% c("Pure"))
dim(merged.roider.singlets)[1] 45068 27981Remove cells with very high nCount_RNA values, set other final QC filters
merged.roider.singlets <- subset(merged.roider.singlets,
subset = nCount_RNA < 40000 & nCount_RNA > 500 & nFeature_RNA > 250)
dim(merged.roider.singlets)[1] 45068 27894Re-compute PCA
Re-scale data now that it has been subset
merged.roider.singlets[["RNA"]] <- split(merged.roider.singlets[["RNA"]], f = merged.roider.singlets$Sample)
merged.roider.singlets <- FindVariableFeatures(merged.roider.singlets,
assay="RNA",
layer="data",
selection.method = "vst",
nfeatures = 5000)
merged.roider.singlets <- ScaleData(merged.roider.singlets, vars.to.regress = "percent.mt")
merged.roider.singlets <- RunPCA(merged.roider.singlets, npcs = 200, verbose = FALSE)Run Harmony integration
integrated.roider <- IntegrateLayers(merged.roider.singlets,
method = HarmonyIntegration,
orig.reduction = "pca",
new.reduction = "integrated.harmony"
)Warning: HarmonyMatrix is deprecated and will be removed in the future from the
API in the futureWarning: Warning: The parameters do_pca and npcs are deprecated. They will be ignored for this function call and please remove parameters do_pca and npcs and pass to harmony cell_embeddings directly.
This warning is displayed once per session.Warning: Warning: The parameter tau is deprecated. It will be ignored for this function call and please remove parameter tau in future function calls. Advanced users can set value of parameter tau by using parameter .options and function harmony_options().
This warning is displayed once per session.Warning: Warning: The parameter block.size is deprecated. It will be ignored for this function call and please remove parameter block.size in future function calls. Advanced users can set value of parameter block.size by using parameter .options and function harmony_options().
This warning is displayed once per session.Warning: Warning: The parameter max.iter.harmony is replaced with parameter max_iter. It will be ignored for this function call and please use parameter max_iter in future function calls.
This warning is displayed once per session.Warning: Warning: The parameter max.iter.cluster is deprecated. It will be ignored for this function call and please remove parameter max.iter.cluster in future function calls. Advanced users can set value of parameter max.iter.cluster by using parameter .options and function harmony_options().
This warning is displayed once per session.Warning: Warning: The parameter epsilon.cluster is deprecated. It will be ignored for this function call and please remove parameter epsilon.cluster in future function calls. Advanced users can set value of parameter epsilon.cluster by using parameter .options and function harmony_options().
This warning is displayed once per session.Warning: Warning: The parameter epsilon.harmony is deprecated. It will be ignored for this function call and please remove parameter epsilon.harmony in future function calls. If users want to control if harmony would stop early or not, use parameter early_stop. Advanced users can set value of parameter epsilon.harmony by using parameter .options and function harmony_options().
This warning is displayed once per session.Transposing data matrixUsing automatic lambda estimationInitializing state using k-means centroids initializationHarmony 1/10Harmony 2/10Harmony converged after 2 iterationsSave integrated object
saveRDS(integrated.roider,"Roider_DLBCL_scRNA_singlets_integrated_harmony.rds")Identify clusters
integrated.roider <- FindNeighbors(integrated.roider, dims = 1:30, reduction = "integrated.harmony")Computing nearest neighbor graphComputing SNNintegrated.roider <- FindClusters(integrated.roider,
resolution = seq(0.05,1,0.05),
algorithm = 2)Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 27948
Number of edges: 873989
Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.9688
Number of communities: 3
Elapsed time: 6 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 27948
Number of edges: 873989
Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.9504
Number of communities: 4
Elapsed time: 6 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 27948
Number of edges: 873989
Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.9368
Number of communities: 6
Elapsed time: 7 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 27948
Number of edges: 873989
Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.9278
Number of communities: 7
Elapsed time: 6 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 27948
Number of edges: 873989
Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.9193
Number of communities: 9
Elapsed time: 7 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 27948
Number of edges: 873989
Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.9113
Number of communities: 9
Elapsed time: 7 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 27948
Number of edges: 873989
Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.9034
Number of communities: 9
Elapsed time: 6 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 27948
Number of edges: 873989
Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.8963
Number of communities: 11
Elapsed time: 6 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 27948
Number of edges: 873989
Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.8895
Number of communities: 12
Elapsed time: 6 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 27948
Number of edges: 873989
Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.8830
Number of communities: 12
Elapsed time: 8 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 27948
Number of edges: 873989
Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.8766
Number of communities: 12
Elapsed time: 6 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 27948
Number of edges: 873989
Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.8701
Number of communities: 12
Elapsed time: 7 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 27948
Number of edges: 873989
Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.8638
Number of communities: 13
Elapsed time: 6 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 27948
Number of edges: 873989
Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.8583
Number of communities: 13
Elapsed time: 6 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 27948
Number of edges: 873989
Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.8531
Number of communities: 14
Elapsed time: 7 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 27948
Number of edges: 873989
Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.8483
Number of communities: 15
Elapsed time: 7 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 27948
Number of edges: 873989
Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.8432
Number of communities: 15
Elapsed time: 6 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 27948
Number of edges: 873989
Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.8384
Number of communities: 16
Elapsed time: 6 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 27948
Number of edges: 873989
Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.8336
Number of communities: 15
Elapsed time: 6 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 27948
Number of edges: 873989
Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.8292
Number of communities: 17
Elapsed time: 7 secondsintegrated.roider <- RunUMAP(integrated.roider,
reduction = "integrated.harmony",
dims = 1:30,
reduction.name = "umap.harmony")11:08:18 UMAP embedding parameters a = 0.9922 b = 1.112Found more than one class "dist" in cache; using the first, from namespace 'BiocGenerics'Also defined by 'spam'11:08:18 Read 27948 rows and found 30 numeric columns11:08:18 Using Annoy for neighbor search, n_neighbors = 30Found more than one class "dist" in cache; using the first, from namespace 'BiocGenerics'Also defined by 'spam'11:08:18 Building Annoy index with metric = cosine, n_trees = 500% 10 20 30 40 50 60 70 80 90 100%[----|----|----|----|----|----|----|----|----|----|**************************************************|
11:08:26 Writing NN index file to temp file /tmp/RtmpGN9GPk/file150e28783fb059
11:08:26 Searching Annoy index using 1 thread, search_k = 3000
11:08:37 Annoy recall = 100%
11:08:39 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
11:08:43 Initializing from normalized Laplacian + noise (using RSpectra)
11:08:43 Commencing optimization for 200 epochs, with 1229528 positive edges
11:09:06 Optimization finishedPlot cluster results across range of resolutions
DimPlot(integrated.roider,
reduction = "umap.harmony",
group.by = paste0("RNA_snn_res.",seq(0.05,1,0.05)),
label = TRUE)
Plot clusters and proceed using res = 0.5 as an example
Idents(integrated.roider) <- integrated.roider$RNA_snn_res.0.5
integrated.roider$seurat_clusters <- integrated.roider$RNA_snn_res.0.5
table(integrated.roider$seurat_clusters)
0 1 2 3 4 5 6 7 8 9 10 11
7348 6288 3436 2410 1906 1900 1869 1085 852 489 278 87 DimPlot(integrated.roider,reduction = "umap.harmony", group.by = "RNA_snn_res.0.5", label = TRUE)
DimPlot(integrated.roider,reduction = "umap.harmony", group.by = "Disease")
DimPlot(integrated.roider,reduction = "umap.harmony", group.by = "Sample")
Plot QC
VlnPlot(integrated.roider, features = "nCount_RNA", group.by="seurat_clusters",pt.size=0) + NoLegend()
VlnPlot(integrated.roider, features = "nFeature_RNA", group.by="seurat_clusters",pt.size=0) + NoLegend()
VlnPlot(integrated.roider, features = "percent.mt", group.by="seurat_clusters",pt.size=0) + NoLegend()
Join with published cluster labels
Note their annotations are provided only on a subset of cells, about 6,000 out of 26,000 total B cells. This means most of our cells have NA for cluster here. We could request cluster labels for all cells from the authors to improve this
overlap <- intersect(rownames(roider.meta$bcells),colnames(integrated.roider))
length(overlap)[1] 5565roider.meta.subset <- roider.meta$bcells[overlap,]
pub.cluster <- roider.meta.subset[colnames(integrated.roider),3]
integrated.roider <- AddMetaData(integrated.roider, pub.cluster, col.name="roider_clusters")
DimPlot(integrated.roider,
group.by = "roider_clusters",
reduction = "umap.harmony",
cells = overlap
) +
DimPlot(integrated.roider,
group.by = "seurat_clusters",
reduction = "umap.harmony",
cells = overlap
)
Plot Roider’s UMAP coordinates for these overlapping cells but colored by our clusters
rownames_to_column(roider.meta.subset,var="Cell") %>%
as_tibble() %>%
inner_join(enframe(integrated.roider$seurat_clusters[overlap]) %>%
dplyr::rename("Cell" = name, "our.cluster" = value),
by = "Cell"
) %>%
ggplot(aes(x=UMAP1,y=UMAP2,color=our.cluster)) +
geom_point(size=0.5) +
theme_classic()
Make Sankey plot comparing cluster labels
roider.clusts <- enframe(integrated.roider$roider_clusters[overlap]) %>%
dplyr::rename("Barcode" = name, "Roider.cluster" = value) %>%
mutate(Roider.cluster = as.numeric(as.character(Roider.cluster)))
our.clusts <- enframe(integrated.roider$seurat_clusters[overlap]) %>%
dplyr::rename("Barcode" = name, "our.cluster" = value) %>%
mutate(our.cluster = as.numeric(as.character(our.cluster)))
make_long (
full_join(roider.clusts,our.clusts,by="Barcode"),
Roider.cluster,our.cluster) %>%
ggplot(aes(x = x,
next_x = next_x,
node = node,
next_node = next_node,
fill = factor(node),
label = node)) +
geom_sankey() +
theme_sankey(base_size = 18) +
theme(legend.position = "none") +
xlab("") +
geom_sankey_text(size = 3, color = "black") +
scale_x_discrete(breaks = c("Roider.cluster","our.cluster"),labels = c("Roider clusters","Our clusters"))
Identify cursory marker genes of each cluster at this resolution
integrated.roider[['RNA']] <- JoinLayers(integrated.roider[['RNA']])
DefaultAssay(integrated.roider) <- "RNA"
vargenes <- presto::wilcoxauc(integrated.roider, 'seurat_clusters', seurat_assay = 'RNA')
top_vargenes <- top_markers(vargenes, n = 100, auc_min = 0.5, pct_in_min = 20, pct_out_max = 20)
top_vargenes# A tibble: 100 × 13
rank `0` `1` `10` `11` `2` `3` `4` `5` `6` `7` `8` `9`
<int> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 1 HSPA… IGLC7 MIR1… GNG11 ATP5… UBE2C CD9 SRGN AC11… GPR1… FBLN5 AL92…
2 2 HSPA… IGHE… EBI3 PCDH9 IGHG3 MKI67 RP5-… LAG3 SELE… TNFR… RP5-… SERP…
3 3 CH17… AL92… IL4I1 GADD… ATP5… NUSA… RP11… SMIM… C1or… LINC… DNTT MAP3…
4 4 RP11… DMD PLEK FBP1 NOP53 CKS1B CCDC… S100… PLAC8 AC24… CCDC… GOLT…
5 5 RP11… MIR3… CD44 CCDC… ATP5… ZWINT YBX3 ATP5… ATP5… CD44 CXor… TCL1B
6 6 TNF ITPR2 NFKB2 ZEB2 SELE… SMC4 KIAA… JSRP1 PLPP5 RIPO… ACY3 CD1C
7 7 RP11… EIF1… MARC… CSNK… SELE… BIRC5 PLAC8 IFI27 AC24… ATP5… CD9 IFNG…
8 8 IGHGP LINC… SLAM… EML6 ELOB TOP2A DNTT MDK AC00… ATP5… KIAA… TYMS
9 9 HIST… BFSP2 BATF BMP3 ATP5… CENPF SOX4 CTSD SLC2… AC11… PLAC8 HTR3A
10 10 SGK1 CCDC… ICAM1 RP11… ATP5… TYMS RASS… LINC… RIPO… PLAC8 XIST SCOC
# ℹ 90 more rowswrite_tsv(top_vargenes,"Roider_DLBCL_scRNA_singlets_integrated_harmony_clustered_prestoMarkers.tsv")Plot top markers identified as a dotplot
top_vargenes <- top_markers(vargenes, n = 5, auc_min = 0.5, pct_in_min = 25, pct_out_max = 10)
top_markers <- top_vargenes %>%
select(-rank) %>%
unclass() %>%
stack() %>%
pull(values) %>%
unique() %>%
.[!is.na(.)]
# Add genes used in B-cell gating above
top_markers <- unique(c(top_markers,"CD79A","CD79B","MS4A1","BANK1","CD3E","CD2","CD3D","GZMK","DERL3", "ITM2C","CTLA4","TIGIT","LAG3","PRF1","SRGN"))
# Compute aggregated expression values of these genes and cluster them to order the figure
rna <- AverageExpression(integrated.roider,assay="RNA",slot="data")As of Seurat v5, we recommend using AggregateExpression to perform pseudo-bulk analysis.
First group.by variable `ident` starts with a number, appending `g` to ensure valid variable names
This message is displayed once per session.rna.sub <- rna$RNA[top_markers,]
cors.genes <- as.dist(1-cor(as.matrix(t(rna.sub)),method="pearson"))
hc.genes <- hclust(cors.genes)
top_markers.sorted <- rownames(rna.sub)[hc.genes$order]
# Plot
DotPlot(integrated.roider,features=top_markers.sorted,assay="RNA",cols=c("blue","red"),cluster.idents=T) + RotatedAxis()
Plot proportions per cluster per disease type
as.data.frame(table(integrated.roider$seurat_clusters, integrated.roider$Sample)) %>%
as_tibble() %>%
dplyr::rename("Sample" = Var2,"Cluster" = Var1) %>%
tidyr::extract(Sample,into=c("Disease","Patient"),regex="(.+)(\\d)",remove=F) %>%
group_by(Sample) %>%
mutate(Proportion = Freq / sum(Freq)) %>%
ggplot(aes(fill = Disease, x = Disease, y=Proportion)) +
geom_boxplot(outlier.shape=NA) +
geom_point(pch = 21, position = position_jitterdodge(jitter.width=0.2)) +
facet_wrap(~Cluster,scales="free")
Compute GCB module score per cell and plot in UMAP space and as dotplot
Derived from HCATonsilData v2 for GCB clusters
combinedGCB <- unique(c("AICDA","MME","RGS13","MEF2B","STMN1","MKI67","PCLAF","NUSAP1","TOP2A","HMGB2","SUGCT","MEF2B","NEIL1","SEC14L1","AICDA","BCL7A","MME","BACH2","TCL1A","CD38","MEF2B","RGS13","SERPINA9","LMO2","HMCES","KLHL6","LRMP","FGD6","MARCKSL1","NEIL1","LMO2","RGS13","MEF2B","BCL2A1","SERPINA9","GMDS","FGD6","RFTN1","PLEK","LRMP","BATF3","MYCBP","NME1","FABP5","PCLAF","MCM7","EBI3","PAICS","PCNA","LMO2"))
integrated.roider <- AddModuleScore(integrated.roider, features = list(combinedGCB), name = "GCB")UMAP
FeaturePlot(integrated.roider, reduction = "umap.harmony", features = "GCB1", order=T) +
FeaturePlot(integrated.roider, reduction = "umap.harmony", features = "HLA-E", order=T) +
plot_layout(ncol=2)
Dotplot
DotPlot(integrated.roider,features=combinedGCB,assay="RNA",cols=c("blue","red"),cluster.idents=T) + RotatedAxis()
Average GCB scores by cluster and sort
as_tibble(rownames_to_column(integrated.roider@meta.data,var="Cells")) %>%
dplyr::select(Cells,seurat_clusters,GCB1) %>%
group_by(seurat_clusters) %>%
summarize(GCB = mean(GCB1)) %>%
dplyr::arrange(desc(GCB))# A tibble: 12 × 2
seurat_clusters GCB
<fct> <dbl>
1 3 0.289
2 4 0.117
3 8 0.103
4 9 0.0605
5 11 0.0451
6 2 0.0201
7 1 0.0102
8 10 -0.00964
9 0 -0.0555
10 6 -0.120
11 5 -0.167
12 7 -0.193 Plot GCB signature vs HLA-E expression
Define which clusters were represented >5% in DLBCL tumors
DLBCLclusters <- as.data.frame(table(integrated.roider$seurat_clusters, integrated.roider$Sample)) %>%
as_tibble() %>%
dplyr::rename("Sample" = Var2,"Cluster" = Var1) %>%
tidyr::extract(Sample,into=c("Disease","Patient"),regex="(.+)(\\d)",remove=F) %>%
group_by(Sample) %>%
mutate(Proportion = Freq / sum(Freq)) %>%
group_by(Disease,Cluster) %>%
summarize(Mean = mean(Proportion)) %>%
dplyr::filter(Disease=="DLBCL") %>%
dplyr::filter(Mean > 0.05) %>%
pull(Cluster) %>%
as.character() %>%
as.numeric()`summarise()` has grouped output by 'Disease'. You can override using the
`.groups` argument.DLBCLclusters[1] 0 1 3 4 5 8Plot boxplot of HLA-E expression divided by GCB vs non-GCB cell clusters, each point is a sample
# Get top scoring cluster for GCB signature
gcb_clust <- as_tibble(rownames_to_column(integrated.roider@meta.data,var="Cells")) %>%
dplyr::select(Cells,seurat_clusters,GCB1) %>%
group_by(seurat_clusters) %>%
summarize(GCB = mean(GCB1)) %>%
slice_max(GCB, n = 1) %>%
pull(seurat_clusters)
as_tibble(rownames_to_column(integrated.roider@meta.data,var="Cells")) %>%
dplyr::select(Cells,Disease,Sample,seurat_clusters,GCB1) %>%
group_by(seurat_clusters,Disease,Sample) %>%
left_join(enframe(integrated.roider@assays$RNA$data["HLA-E",],name = "Cells",value = "HLAE") %>% as_tibble(), by="Cells") %>%
summarize(GCB = mean(GCB1),HLAE = mean(HLAE)) %>%
dplyr::filter(Disease=="DLBCL" & seurat_clusters %in% DLBCLclusters) %>%
mutate(CellClass = case_when(
seurat_clusters == gcb_clust ~ "GCB",
T ~ "NonGCB"
)
) %>%
ggplot(aes(x=CellClass,y=HLAE,fill=CellClass)) +
geom_boxplot(outlier.shape=NA) +
geom_point(pch = 21, position = position_jitterdodge(jitter.width=0.2))`summarise()` has grouped output by 'seurat_clusters', 'Disease'. You can
override using the `.groups` argument.
Repeat above but plot GCB from rLN vs all other cells
as_tibble(rownames_to_column(integrated.roider@meta.data,var="Cells")) %>%
dplyr::select(Cells,Disease,Sample,seurat_clusters,GCB1) %>%
left_join(enframe(integrated.roider@assays$RNA$data["HLA-E",],name = "Cells",value = "HLAE") %>% as_tibble(), by="Cells") %>%
mutate(CellClass = case_when(
seurat_clusters == gcb_clust & Disease == "rLN" ~ "GCB",
seurat_clusters == gcb_clust & Disease != "rLN" ~ "GCBother",
T ~ "NonGCB"
)
) %>%
group_by(CellClass,Disease,Sample) %>%
summarize(GCB = mean(GCB1),HLAE = mean(HLAE)) %>%
ggplot(aes(x=Disease,y=HLAE,fill=Disease)) +
geom_boxplot(outlier.shape=NA) +
geom_point(pch = 21, position = position_jitterdodge(jitter.width=0.2)) +
geom_text(aes(label=Sample),hjust=0.5,vjust=0.5) +
facet_grid(~CellClass,scales="free_x",space = "free_x")`summarise()` has grouped output by 'CellClass', 'Disease'. You can override
using the `.groups` argument.
Plot scatterplot of GCB signature vs HLA-E expression, colored by GCB vs non-GCB, each point is a cell
as_tibble(rownames_to_column(integrated.roider@meta.data,var="Cells")) %>%
dplyr::select(Cells,Disease,Sample,seurat_clusters,GCB1) %>%
group_by(seurat_clusters,Disease,Sample) %>%
left_join(enframe(integrated.roider@assays$RNA$data["HLA-E",],name = "Cells",value = "HLAE") %>% as_tibble(), by="Cells") %>%
dplyr::filter(Disease=="DLBCL" & seurat_clusters %in% DLBCLclusters) %>%
mutate(CellClass = case_when(
seurat_clusters == gcb_clust ~ "GCB",
T ~ "NonGCB"
)
) %>%
ggplot(aes(x = GCB1, y = HLAE, color = CellClass)) +
geom_point(cex=0.5) +
geom_smooth(method = "lm") +
theme_classic() +
xlab("GCB signature score") +
ylab("HLA-E expression") +
scale_color_manual(values=c("#B07AA1","#499894")) +
stat_cor(aes(color = CellClass), method = "pearson", label.x = 0.2)`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'png
2 Plot ridgeline showing HLA-E expression divided by GCB vs non-GCB
as_tibble(rownames_to_column(integrated.roider@meta.data,var="Cells")) %>%
dplyr::select(Cells,Disease,Sample,seurat_clusters,GCB1) %>%
group_by(seurat_clusters,Disease,Sample) %>%
left_join(enframe(integrated.roider@assays$RNA$data["HLA-E",],name = "Cells",value = "HLAE") %>% as_tibble(), by="Cells") %>%
dplyr::filter(Disease=="DLBCL" & seurat_clusters %in% DLBCLclusters) %>%
mutate(CellClass = case_when(
seurat_clusters == gcb_clust ~ "GCB",
T ~ "NonGCB"
)
) %>%
ggplot(aes(x = HLAE, y = CellClass, fill = CellClass)) +
geom_density_ridges(alpha=0.7) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_discrete(expand = c(0, 0)) +
coord_cartesian(clip = "off") +
theme_ridges(grid = FALSE, center_axis_labels = TRUE)Picking joint bandwidth of 0.127
Save clustered object
saveRDS(integrated.roider,"Roider_DLBCL_scRNA_singlets_integrated_harmony_clustered.rds")Get session info
sessionInfo()R version 4.3.1 (2023-06-16)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Rocky Linux 8.10 (Green Obsidian)
Matrix products: default
BLAS/LAPACK: /usr/lib64/libopenblasp-r0.3.15.so; LAPACK version 3.9.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
time zone: America/New_York
tzcode source: system (glibc)
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] ggpubr_0.6.0 ggsankey_0.0.99999
[3] scGate_1.6.0 ggridges_0.5.6
[5] scDblFinder_1.14.0 harmony_1.2.0
[7] vctrs_0.6.5 patchwork_1.3.0
[9] scater_1.34.0 scuttle_1.10.3
[11] speckle_1.0.0 Matrix_1.6-4
[13] readxl_1.4.3 SingleCellExperiment_1.22.0
[15] SummarizedExperiment_1.30.2 Biobase_2.60.0
[17] GenomicRanges_1.52.1 GenomeInfoDb_1.36.4
[19] IRanges_2.34.1 S4Vectors_0.38.2
[21] BiocGenerics_0.46.0 MatrixGenerics_1.12.3
[23] matrixStats_1.2.0 flexmix_2.3-19
[25] lattice_0.22-5 SeuratWrappers_0.3.19
[27] miQC_1.8.0 lubridate_1.9.3
[29] forcats_1.0.0 stringr_1.5.1
[31] dplyr_1.1.4 purrr_1.0.2
[33] readr_2.1.5 tidyr_1.3.1
[35] tibble_3.2.1 ggplot2_3.5.1
[37] tidyverse_2.0.0 Seurat_5.1.0
[39] SeuratObject_5.0.2 sp_2.1-3
[41] sctransform_0.4.1 glmGamPoi_1.12.2
[43] presto_1.0.0 Rcpp_1.0.12
[45] devtools_2.4.5 usethis_2.2.2
[47] data.table_1.15.0
loaded via a namespace (and not attached):
[1] fs_1.6.3 spatstat.sparse_3.0-3
[3] bitops_1.0-7 httr_1.4.7
[5] RColorBrewer_1.1-3 backports_1.4.1
[7] profvis_0.3.8 tools_4.3.1
[9] utf8_1.2.4 R6_2.5.1
[11] mgcv_1.9-1 lazyeval_0.2.2
[13] uwot_0.1.16 urlchecker_1.0.1
[15] withr_3.0.0 gridExtra_2.3
[17] progressr_0.14.0 cli_3.6.2
[19] spatstat.explore_3.2-6 fastDummies_1.7.3
[21] labeling_0.4.3 spatstat.data_3.0-4
[23] pbapply_1.7-2 Rsamtools_2.16.0
[25] R.utils_2.12.3 parallelly_1.37.0
[27] sessioninfo_1.2.2 limma_3.56.2
[29] generics_0.1.3 BiocIO_1.10.0
[31] vroom_1.6.5 ica_1.0-3
[33] spatstat.random_3.2-2 car_3.1-2
[35] ggbeeswarm_0.7.2 fansi_1.0.6
[37] abind_1.4-5 R.methodsS3_1.8.2
[39] lifecycle_1.0.4 yaml_2.3.8
[41] edgeR_3.42.4 carData_3.0-5
[43] SparseArray_1.2.2 Rtsne_0.17
[45] grid_4.3.1 promises_1.2.1
[47] dqrng_0.3.2 crayon_1.5.2
[49] miniUI_0.1.1.1 beachmat_2.16.0
[51] cowplot_1.1.3 metapod_1.8.0
[53] pillar_1.9.0 knitr_1.45
[55] rjson_0.2.21 xgboost_1.7.7.1
[57] future.apply_1.11.1 codetools_0.2-19
[59] leiden_0.4.3.1 glue_1.7.0
[61] remotes_2.4.2.1 png_0.1-8
[63] spam_2.10-0 cellranger_1.1.0
[65] gtable_0.3.4 cachem_1.0.8
[67] xfun_0.42 S4Arrays_1.2.0
[69] mime_0.12 survival_3.5-8
[71] statmod_1.5.0 bluster_1.10.0
[73] ellipsis_0.3.2 fitdistrplus_1.1-11
[75] ROCR_1.0-11 nlme_3.1-164
[77] bit64_4.0.5 RcppAnnoy_0.0.22
[79] irlba_2.3.5.1 vipor_0.4.7
[81] KernSmooth_2.23-22 colorspace_2.1-0
[83] nnet_7.3-19 ggrastr_1.0.2
[85] UCell_2.6.2 tidyselect_1.2.0
[87] bit_4.0.5 compiler_4.3.1
[89] BiocNeighbors_1.18.0 DelayedArray_0.26.7
[91] plotly_4.10.4 rtracklayer_1.60.1
[93] scales_1.3.0 lmtest_0.9-40
[95] digest_0.6.34 goftest_1.2-3
[97] spatstat.utils_3.0-4 rmarkdown_2.25
[99] XVector_0.40.0 htmltools_0.5.7
[101] pkgconfig_2.0.3 sparseMatrixStats_1.12.2
[103] fastmap_1.1.1 rlang_1.1.3
[105] htmlwidgets_1.6.4 shiny_1.8.0
[107] DelayedMatrixStats_1.22.6 farver_2.1.1
[109] zoo_1.8-12 jsonlite_1.8.8
[111] BiocParallel_1.34.2 R.oo_1.26.0
[113] BiocSingular_1.16.0 RCurl_1.98-1.14
[115] magrittr_2.0.3 modeltools_0.2-23
[117] GenomeInfoDbData_1.2.10 dotCall64_1.1-1
[119] munsell_0.5.0 viridis_0.6.5
[121] reticulate_1.35.0 stringi_1.8.3
[123] zlibbioc_1.46.0 MASS_7.3-60.0.1
[125] plyr_1.8.9 pkgbuild_1.4.3
[127] parallel_4.3.1 listenv_0.9.1
[129] ggrepel_0.9.5 deldir_2.0-2
[131] Biostrings_2.68.1 splines_4.3.1
[133] tensor_1.5 hms_1.1.3
[135] locfit_1.5-9.8 igraph_2.0.2
[137] spatstat.geom_3.2-8 ggsignif_0.6.4
[139] RcppHNSW_0.6.0 reshape2_1.4.4
[141] ScaledMatrix_1.8.1 pkgload_1.3.4
[143] XML_3.99-0.16.1 evaluate_0.23
[145] scran_1.28.2 BiocManager_1.30.22
[147] tzdb_0.4.0 httpuv_1.6.14
[149] RANN_2.6.1 polyclip_1.10-6
[151] future_1.33.1 scattermore_1.2
[153] rsvd_1.0.5 broom_1.0.5
[155] xtable_1.8-4 restfulr_0.0.15
[157] RSpectra_0.16-1 rstatix_0.7.2
[159] later_1.3.2 viridisLite_0.4.2
[161] memoise_2.0.1 beeswarm_0.4.0
[163] GenomicAlignments_1.36.0 cluster_2.1.6
[165] timechange_0.3.0 globals_0.16.2